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^ ' Abstract 

in : 

b^ . We implement fermions on dynamical random triangulation 

O '■ and determine numerically the spectrum of the Dirac- Wilson op- 

^ I erator V for the system of Majorana fermions coupled to two- 

■^ I dimensional Euclidean quantum gravity. We study the depen- 

^ ' dence of the spectrum of the operator eV on the hopping param- 

A . eter. We find that the distributions of the lowest eigenvalues be- 

D ! come discrete when the hopping parameter approaches the value 

"T^ I l/v3. We show that this phenomenon is related to the behavior 

.^ ' of the system in the 'antiferromagnetic' phase of the correspond- 

^ ' ing Ising model. Using finite size analysis we determine critical 

j^ . exponents controlling the scaling of the lowest eigenvalue of the 

spectrum including the Hausdorff dimension dn and the exponent 
K which tells us how fast the pseudo-critical value of the hopping 
parameter approaches its infinite volume limit. 

Introduction 

The dynamical triangulation approach to quantum gravity has proven to be 
a very powerful method fl], ||, ^. In two- dimensions it yields the same results 
for critical exponents as the Liouville theory 0, ^, ^ . Contrary to the latter, 
this approach can be straightforwardly generalized to higher dimensional 
case which is frequently referred to as simplicial gravity 0, H. Results from 
numerical studies of pure gravity without matter fields in four dimensions 
showed that the continuum limit of this model does not exist iPl. In order 



to obtain more realistic models, one has tried to include matter fields and 
to couple them to gravity ||10[- This program has so far succeeded only for 



bosonic matter. Putting fermions on random simplicial manifold is a more 
difficult task. In general it requires introducing an additional field of local 



frames in order to define a spin structure |Tl], |12|, ^. In the case of a 
compact manifold this is a topological problem. Although many ingredients 
of the construction are known and can be generalized to any number of 
dimensions, the topological part of the problem has been solved so far only 
in two dimensions [|13], Q. 

In this paper we will study properties of the Dirac- Wilson operator on 
two-dimensional dynamical triangulation with spherical topology. The anal- 
ysis of the spectrum in the critical region allows us to calculate critical indices 
as for example the Hausdorff dimension. 

We cross-check properties of the spectrum using the fact that the partition 
function of the fermionic model can be mapped into the partition function 
of Ising model on dynamical triangulation, which is analytically solvable 



151, |T6l|T7i 



The spectrum of the Major ana- Dirac- Wilson operator eD becomes dis- 
crete when the hopping parameter admits the value l/-\/3 corresponding to 
the value /3 = of the coupling in the Ising model. We show that this 
behaviour can be explained by the presence of a set of points in 'antifer- 
romagnetic' phase (/3 < 0), for which some eigenvalues of the operator are 
determined by local properties of the triangulation and not by its random 
character. 

The paper is organized as follows : First we define the model, then we 
recall some facts about its relation to the Ising model [0, we present results 
of numerical studies and shortly conclude at the end by summarizing and 
listing open questions. In the appendix, for comparison, we calculate the 
spectrum of the Dirac- Wilson operator on a regular triangulation. 

The model 

The model of fermions minimally coupled to Euclidean gravity is given by 
the partition function 






TeT TeT 



where the sum goes over d-dimensional simplicial manifolds from a class T, 
say, for instance, with spherical topology. Each triangulation is dressed with 
the fermion field located in the centers of d-simplices. The integral over field 



on a given triangulation T defines the partition function Zt^ wliicli at tfie 
same time provides a weigfit of tfiis triangulation in tlie ensemble. The action 
reads 

St = -KY, ^i^.,^, ^\Y, '^^^^ = Y. ^^^^^-^^^ ' (2) 

where the fermionic fields \E'j are located in the centers of triangles. The sum 
over (zj) runs over oriented pairs of neighboring triangles, or equivalently, 
over oriented dual links. The hopping operator is 

^.. = ^(l + 4^-7)W.,. (3) 

The Dirac- Wilson operator is denoted by T>ij and the spin connection by Uij. 
In order to be able to calculate spinor and vector components, we endow each 
(i-simplex with an orthonormal local frame. A frame is a set of orthonormal 
oriented vectors Cq, a = 1, . . . c?. To each vector Ca we ascribe a Dirac gamma 
matrix 7", in such a way that its numerical value is identical in each frame. 
The local vector iiij in eq. (^) is a unit vector which points from the center 
of the simplex j to the center of one of its neighbors i. It just tells us the 
direction of the local derivative. The inner product of this vector and of 
gamma matrices, which is denoted by dot in (0), has to be understood as a 
sum of gamma matrices 7" multiplied by the components of \fii}\a in the given 
frame at z, denoted by the upper index. Thus, the product of the same vector 
riij expressed in another frame yields a different matrix : jt-^^ ■ 7 = ['t-j-'^ ]a7"- 

As mentioned the matrix Uij plays the role of spin connection. It allows 
us to parallel transport a spinor from the simplex j to the simplex i, or 
in other words, to recalculate spinor components between two neighboring 
frames i and j. The matrix lAij is an image in the spinorial representation of 
the rotation matrix Vij which parallel transports vectors. The map Vij -^ Uij 
is not unique, namely it is defined up to sign. As we will see below, the signs 
oiU must be adjusted to fulfill a consistency condition (H) for all elementary 
plaquettes of the simplicial manifold. This is a topological problem. 

This problem has been solved in two-dimensions where an explicit con- 
struction of the signs of the spin connection matrices Uij has been given |jl3l . 
Let us shortly recall the main steps of the construction. 

In two dimensions each orthonormal frame consists of two vectors Cja 
where a is 1 or 2. The first index of e^a refers to the triangle in which the 
frame is located. For any pair of neighboring triangles i, j we can define 
a spin connection as a two by two rotation matrix [t/jj]^, such thatQ tia = 

^In general a connection can be a dynamical field. 



Tl,bWii\^a^ib- Using matrix notation this relation can be written as Cj = t/jjCj, 
where 

U.. = ^^^<t>^, = ( cosA0,j- sinA0ij- \ 
'■^ \^ — sin A0jj cos A0jj J 

and A0jj is the relative angle between the two neighboring frames, e is the 
standard antisymmetric tensor. 

The trace of an elementary loop around a dual plaquette is a geometrical 
invariant directly related to the curvature (deficit angle) of the vertex in the 
center of the plaquette. One can check that 

-TiUU ...U = -Tr e^(2^-^^) = cos Ap , (5) 

where Ap is the deficit angle of the vertex in middle of the plaquette. The 
product UU . . .U oi connections on all links on the plaquette perimeter P 
is a rotation matrix which gives the integrated rotation of a tangent vector 
parallel transported around this loop. The equation (|) is a sort of Wil- 
son discretization |jT8|, |19[ of curvature calculated from the Cartan structure 
equations [^ . 

Now the idea is to write down an analogous equation as (|^) in the spino- 
rial representation. First we have to introduce a parallel transporter Uij for 
spinors for each pair of neighboring vertices. This is exactly the object which 
we need in (j^). The connection Uij is an spinorial image of Uij. One can 
choose a representation of gamma matrices such that Uij = U^-. One immedi- 
ately sees that indeed Uij can be calculated for a given Uij up to sign. When 
defining the Dirac- Wilson operator (^ we cannot allow for ambiguities, so 
we have to give a unique prescription how to calculate Uj. We do this by 
choosing 

...,, / -s^ sin^\ 
U,, = e-^=\ (6) 



_sin^ cos^ 



and specifying the angles A0jj uniquely. More precisely we define A0jj = 
Vi ~ yj + ^ where the angle cpf at triangle j is the angle between the 
vector Cji of the frame at j and the vector riij (pointing from j to i), and 
likewise cpj at triangle i is the angle between the vector en and the vector riji 
(pointing from i to j) (see fig.|l|). Both the angles are restricted to the range 
[0, 2tt) and both are measured in the same direction, say clockwise. Thus 
the angle A0jj is defined without the 2tt ambiguity and hence the rotation 
matrix U^j is also uniquely determined including the total sign. 




n 



Figure 1: Local geometry of two neighboring triangles is shown. The po- 
sition of the first frame vector ei for a given triangle is marked by a line 
emerging from the triangle center. The position of the second frame vector 
62 is implicitly given by the fact that the angle between ei and 62 counted 
clockwise is 7r/2. The vector riji points from the center of the triangle i to 
j. The arch in the triangle i represents the angle (pj between en and riji. 
The arch in the triangle j corresponds to the angle 0^- between Cji and Uij. 



In the example shown in figure 
neighbors of j : (f))^ = tc/3, tpn - 



4" 

7r. 



vr, 0^ 
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Svr/S, and for other two 



One can easily check that, for the definition (^) of Wy's, the parallel 
transporter around an elementary loop gives 



-TiUlA . . M = Sp cos 



A, 



(7) 



The argument of cosine got halved Ap -^ Ap/2 in comparison with (^) 
because for each link on P we have U^ = U . The total sign Sp of the 
product UU . . .U has to be calculated. It turns out to depend on all angles 
A0jj on the loop and it may admit either value ±1 |[T3| . 

The presence of elementary plaquettes which would have negative sign 
is an unwanted effect. For example, a spinor transported around a flat pla- 
quette with Sp = —1 would change the sign ip -^ —ip. We require that the 
parallel transport around a close loop in a fiat patch not change a spinor. 
Furthermore, we require the signs Sp to be positive for all elementary pla- 
quettes 



S, 



VP. 



One can give the following argument in favor of the naturalness of this re- 
quirement. An elementary loop goes through triangles sharing a vertex. The 
geometry of a patch consisting of those triangles corresponds to the geometry 
of a cone. It is everywhere fiat except at the vertex where it is singular. One 
can regularize such geometry by smoothing the peak of the cone (making it 
differentiable) in a very small region with radius e ~ 0. Such a regulariza- 
tion does not affect the loop which lies in a distance R^ e from the vertex. 
Continuously shrinking the loop in such a regularized geometry, one can con- 
tinuously change the angle of the loop rotation matrix without changing the 
sign. A completely shrunken loop must have positive sign since it lies in a 
flat patch. This implies Sp = +1. One can also check that the consistency 
condition (^) plays an essential role in the topological considerations or in 
deriving the equivalence with the Ising model. 

The construction of the connections given in (^) does not fulflll the con- 
sistency condition (^. We will therefore modify the construction of Ws by 
introducing for each link an additional sign degree of freedom Sij 

Uij = SijC 2 . (9) 

One can show that this freedom is sufficient to globally, for each elementary 
loop, fulflll the consistency condition (H), on a triangulation of an orientable 
manifold. Thus, technically, to deflne the Dirac- Wilson operator on a trian- 
gulation, we have to flrst assign an orthonormal frame to each triangle, and 
then for the frame assignment, to flnd link signs Sij meeting the consistency 
condition (|]) for each plaquette. The remaining part is straightforward. 
Namely, we express the operator U^j in terms of the angles (pl and 0^ and 
likewise, the product nlj ■ 7 in terms of 0^ . Thus, we parameterize the 

hopping operator Hij entirely by 0^ and 0^ and Sij. For each pair of neigh- 
boring triangles the angles can be read off from the given frame assignment 
(see flg.|l]). 

Choosing the following representation of gamma matrices : 71 = cts, 
72 = CTij where the a are Pauli matrices, we eventually arrive at 



sm ^- cos -^ sin -^ sin -^ 
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(10) 



We see that in two dimensions the Dirac- Wilson operator on a triangulation 
T is given by a matrix consisting of two-by-two blocks 



where Pij is an adjacency matrix : 

1, if i,j are neighbors 



^■' I 0, otherwise 

There is no summation over ij in the equation (0). The blocks Tiij have a 
very simple structure. In fact, we can simplify it further by restricting the 
set of values of the angles in (|10|) from the whole interval [0, 27r) to a discrete 
set of three values separated by 27r/3, for instance, 7r/3,7r, Stt/S. For this 
choice, the first vector ei of a frame at a triangle points from the center of 
the triangle to one of its vertices. In a sense, this set of three frame positions 
is a minimal set reflecting the symmetry of equilateral triangle. 

Since physical quantities do not depend on the choice of frames, this 
restriction is a sort of gauge condition. With this choice, the blocks ([10| ) may 
admit only nine different forms depending on the nine different choices of the 
angles 0^, 0* in (|10D . They can be precomputed. For example, for the frame 

assignment in fig.[5 (pf = Stt/S, (pj = tx and for Sij = 1 we have 



^ i \ .. / 
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Fermions and the Ising model 

The idea is now to calculate spectra of the Dirac- Wilson operator for different 
triangulations from the ensemble (|I|). Summing up (averaging) all the spectra 
we obtain the spectrum of the Dirac- Wilson operator for fermions interacting 
with 2d gravity. More precisely, we will consider a field of Majorana-fermions 
coupled to gravity. At the critical point it corresponds to the conformal field 
with the central charge c = 1/2. 

Denote the components of the spinor \1/ by \I/q,, and of ^ by \I'^. The 
Majorana condition reads : \1/^ = e^°"^a or ^a = '^^^/3a, where e is the 
standard antisymmetric tensor. In this notation, the action for Majorana 
fermions can be written as 

st = J2 *r [%]^^./^ = E ^-^:^'^.73 ' (14) 



y 



where 



Vff = e-^[P.,]? = ^e'^%, - irP,,e-^[Hd? , (15) 



or in short T> = eD : 

V = Vo-Kn (16) 

where Pq = |e x 1 is a deterministic part, and 7i is a random part which 
consists of two-by-two matrices Tiij = eTiij (^U|) located at nonvanishing 
positions of the adjacency matrix Pij which correspond to pairs of neighboring 
triangles. Pij is an off-diagonal random matrix in the ij indices. 

One can show that V is antisymmetric under the change of pairs of indices 

^tf = -^J" , (17) 

and hence Pf P = DetP = Detl^. For each triangulation individually, 
the integral over fermions in (|I]) yields Pfaffian of the matrix V. Thus for 
Majorana fermions on a two-dimensional triangulation the partition function 
(|l|) is a sum of Pfaffians of the Dirac- Wilson operator 

In the last step we have used the inequality PiV^ > 0, which can be proven 
by the hopping parameter expansion. The consistency condition (§) turns 
out to be essential in the proof. Namely, one shows that the PfFafian is 
represented as a sum over loop configurations each of which contributes a 
positive factor if the condition (|) is met |TB|, |n|. 



Using this expansion one can also establish the equivalence between the 
partition function Zj- and the partition function of the nearest neighbor Ising 
model with spins o"^ located at the vertices 2=,, of T 

Zt= J2 ""'^^^ (^9) 

W*}t 



where 



Et = -J2 K^>-1)- (20) 

The partition functions Z^ (|T9D and Z^ (|l]) are equal for 

In the derivation of this equivalence one identifies loop configurations, arising 
in the hopping expansion of Zt with domain walls of the Ising model. Again, 

8 



the consistency condition @ plays the crucial role here. For a non-spherical 
triangulation one has to carefully treat topological effects related to the ex- 
istence of non-contractable loops which may give a negative contribution for 
antiperiodic boundary conditions. One can get rid of all negative contribu- 
tions performing the GSO projection that is summing over all spin structures 
of the manifold [^]. This and another topological issues are discussed else- 



where flj]. Here we will restrict ourselves to spherical triangulations for 
which we automatically have Zj- = Z^ for each triangulation T ^ T and 
hence also for the sum over all triangulations in T 

Z = Z = ^Zt. (22) 

tgt 

The critical temperature of the Ising model for the partition function Z is 



known analytically to be /3 = ^ In ^ ^ 0.2162730 |jT^. Translating the 



critical temperature to the hopping parameter (^) we obtain the following 
critical value 

K^^ = ^^ ^ 0.3746166 (23) 

393 ^ ^ 

for which fermions become massless. Another interesting point which can be 
deduced from the equation (^Tj) is that the value of the hopping parameter 
Kq = l/y/S corresponds to /3 = which is the border between the ferromag- 
netic and antiferromagnetic regimes. For /3 < one expects frustration in 
the Ising model on a triangulation and hence that the lowest energy state is 
highly degenerated. As we will see below, the behavior of the spectrum of 
the Dirac- Wilson is also sensitive to passing over this border. First, we will 
however restrict K to the 'ferromagnetic' range [0,-K'o]. 

The equivalence of the partition functions Zt and Zt may be used to 
relate the average energy Et of the Ising model on a triangulation T to 



eigenvalues of the Dirac- Wilson operator. Differentiating both sides of (|T8 
with respect to (3 we obtain 

where Aq are eigenvalues of the Dirac- Wilson operator Vt- For our choice 
of the representation of gamma matrices, Vt is a real matrix. Its spectrum 
consists of either real eigenvalues or of pairs of complex conjugates. Thus 
the sum on the right hand side of (p^ is a real number. Similarly, the 
fluctuations of the Ising energy on the triangulation T are given by 

^^- A^ -nW^~ 2A^ + ^^^- ^^^^ 
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Figure 2: Heat capacity c„(/5) for the system with A^ = 16 triangles for the 
three cases discussed in the text : (a) ising-ising (hne); (b) ising-fermion 
(filled symbols); (c) fermion-fermion (empty symbols). The three methods 
give the same results within the error bars, which are here of order of the 
size of the symbols used. 



Averaging over triangulations we obtain the energy density and heat-capacity 
of the Ising model coupled to gravity calculated in terms of the eigenvalues 
of the Dirac- Wilson operator 



{er) = 1 



2N ^ " 



(26) 



/32(4) = P'' 



[-{^EK^)^{^EK')] 



(27) 



The equivalence of the models can also be very useful in MC simulations 
of the model. To show this, let us compare three numerical experiments in 
which (a) the Ising model is used to generate triangulations and to measure 
the Ising energy and heat capacity; (b) the Ising model is used as a generator 
of triangulations but measurements are carried out using the fermion field; (c) 
the fermionic model is used both to generate triangulations and to perform 
measurements. As shown in fig.0 the three methods yield exactly the same 
results. The methods differ, however, significantly in the CPU time needed 
to generate results of the same quality. The first difference comes from the 
configuration generator which is much faster for the Ising model than for the 
fermionic determinant. In the latter case, to calculate a Metropolis weight 



10 



for a single local change of triangulation, i.e. a flip of one link on the trian- 
gulation, requires the recomputation of the determinant of the Dirac- Wilson 
operator on the modified lattice. This is a tedious task for which the number 
of operations grows with the third power of the system size N. Thus one 
expects that the time of a sweep through the lattice grows as N"^ for the 
fermionic configurations generator. One sweep for the Ising model, which 
consists of a sweep of local updates of Ising spins, a fixed number of Wolff 
cluster updates, and a sweep of local changes of triangulation, lasts in CPU 
units roughly proportionally to the system size A^. Thus, the fermionic al- 
gorithm is competitive with the Ising generator only for very small lattices. 
As far as measurements are concerned the situation is more complex. For 
example, one cannot determine the spectrum of the Dirac- Wilson operator 
using only the Ising spins. One can, however, do the opposite. For a given 
lattice, the time of calculating all eigenvalues of the Dirac- Wilson operator 
is proportional to A^^. Having done this, one is able for this triangulation to 
exactly calculate the Ising energy (|2^ ) and it higher moments (pSD without 
statistical fluctuations. If one instead used the Ising model, one has to sam- 
ple Ising configuration many times to reduce the error. In general, the cost 
of a single measurement of the energy is proportional to A^. The error of the 
single measurement of the energy density decreases like I/a/ZV. Summariz- 
ing, we expect the CPU time to measure energy with a given precision to 
grow as y/N. The CPU time grows rapidly with the order for measurements 
of higher moments of energy. 

In order to obtain the data the quality presented in fig. |] for A^ = 16, 
the methods discussed above required (a) 1000 CPU min. (b) 6 CPU min. 
and (c) 100 min. on the computer Alpha XP1000/EV6/500 MHz. 

Spectrum of the Dirac- Wilson operator 

In the production runs we use the method (b), which relies on generating 
triangulations from the partition function of the Ising model. At each mea- 
surement we ignore the values of the Ising spins and we assign frames Cj 
and Sjj-signs to the triangulation to reconstruct the Dirac- Wilson operator 

(00) • 

A typical spectrum of the Dirac- Wilson operator on random triangulation 

is shown in fig.^. The main effect on the spectrum of changing the hopping 

parameter K is to rescale it around the point (|,0). The positions of the 

two claw-shaped ends of the spectrum move with K. One can find a value of 

K for which the ends lie closest the origin (0,0). This value can be treated 

as a pseudo-critical value K^, for which the mass of the fermion excitation 
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Figure 3: The distribution of eigenvalues A of the Dirac- Wilson operator for 
A^ = 64 and K = 0.364 on random lattice. 



is minimal. For K < K^, the origin (0, 0) lies outside the claws, while for 
K > K^: inside. In fact, this is the main difference between the two regimes, 
since beside the scaling factor the shape of the spectrum is almost constant. 

The claw-shaped ends of the pseudo-critical spectra successively approach 
each other when the system size A^ is increased. They eventually close en- 
tirely at the origin (0, 0) for infinite A^, signaling the occurrence of massless 
excitations. 

One can alternatively study the spectrum of the operator P = eP. In 
fact, this operator is closer related, in spirit, to Majorana fermions (|15|) due 
to the presence of the charge conjugation matrix e. Since the matrix V is 
antisymmetric, and it is real in our representation, its spectrum is purely 
imaginary. Thus, the eigenvalue density of this operator is one dimensional : 



p(x) 



hm -- 



Y^6{x-zX)\ 



{2i 



The spectrum is symmetric p(x) = p(— x). We will constrain ourselves to 
the positive branch. For each triangulation eigenvalues of the positive part 
of the spectrum can be ordered : Aq < Ai < A2 < . . . . Collecting separately 
histograms for the lowest, second lowest, third lowest eigenvalues e.t.c. one 
obtains distributions Pj(x) of the j-th eigenvalue. Of course, by construction : 

We studied numerically the dependence of pj (x) on the hopping param- 
eter. In figure ^ is shown the distribution ^o{x) for different values of K. 
One can see that it is discrete for Kq = l/\/3 consisting of separate nar- 
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Figure 4: Evolution of the shape of the probabihty distribution 'po{x) of the 
lowest eigenvalue of the operator eV for A^ = 64 as a function of the hopping 
parameter K : {&) K = 1/^3 = 0.5774, (b) K = 1/V3e~°-^ = 0.4727, (c) 
K = l/V3e-°-^ = 0.3870, (d) K = 1/V^e-^-^ = 0.3169. Each histogram 
presented in the figure contains 7 x 10^ counts. The bin size is 5 x 10"^. The 
histograms are normalized to have area one. 



row peaks. This may appear surprising at first glance. However, when K 
becomes smaller than Kq this spectrum becomes continuous : it gradually 
changes when A = Kq — K grows to eventually become a smooth bell-shaped 
distribution (fig.§). As we shall see below the discreteness of the spectrum 
at Kq = l/v3 is not a finite size effect. We made the following experiment. 
We generated a quenched ensemble of random triangulations by ignoring the 
coupling of fermions to gravity. It is an ensemble of triangulations for pure 
gravity. Then, for each triangulation from this ensemble, we calculated the 
operator (|15D and determined its lowest eigenvalues for different values of K. 
One can expect that outside the critical region, where fermions are massive, 
the approximation should not significantly affect the spectrum of the model. 
Indeed we checked numerically for a few values of K that it is practically 
impossible to distinguish between the spectrum for the quenched and the 
full model. 

The results of the quenched experiment are presented in figure ^ where 
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Figure 5: (a) Evolution of the lowest eigenvalue of the operator (p!5|) with 
the hopping term, calculated on an ensemble of given triangulations with 
A^ = 32 triangles, (b) Evolution of seven lowest eigenvalues of the operator 
([T3|). The cross-points of the bundles in the 'antiferromagnetic' phase are 
seen. The cross-points can be numbered by successive integers from 6 to 13 
corresponding to the length q of the elementary loops, as discussed in the 
text. 



one can see lowest eigenvalues of the operator (|T5|) for different K. The 
'ferromagnetic' region corresponds to the interval < i^ < 1/a/3 = 0.5774, 
and the 'antiferromagnetic' one K > 1/^/3. Universal properties of the 
model are encoded in the behavior of the spectrum around the critical value 
Kcr = 0.3746 (^), lying deep inside the 'ferromagnetic' phase, where the 
eigenvalue bundle has a dip. For K -^ only the deterministic part |e"^(5jj 
of the operator (|15D survives. 

The bundles of eigenvalues have an interesting property. The lines of the 
bundle cross at some points in the 'antiferromagnetic' phase (see fig. U). The 
meaning of a cross-point is that, for the corresponding value of K, the opera- 
tor (|15]) has a common identical eigenvalue for many different triangulations. 
The reason why it happens is that this eigenvalue is entirely related to the 
existence of elementary loops of the length q on the triangulation and not 
to the whole random structure of the triangulation, as it generically takes 
place. 

In order to understand the mechanism of the occurrence of the cross- 
points, consider an elementary loop on the dual lattice of length q consisting 
of vertices ii, i2, ■ ■ ■ , iq- One can show that the corresponding submatrix 2g x 
2g of ( |T5| ) built of the two-by-two blocks at the (zi, 22, • • • , iq) x (^1, "^2, • • • , ig) 
positions has for some Kg an eigenvalue Xg which depends only on q. The 
pair {Kg, Xg) corresponds to a cross-point in the figure (|^). Furthermore, it 



turns out that for K 



Kg and A 



Xg the entire 2g rows for ii,i2, 
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of the matrix I? — A^l are linearly dependent. This means that Xg is not 
an eigenvalue of the submatrix but also of the whole matrix V for K = 



Kq (|T5|). The pairs {Kg, Xq) correspond to the cross-points of eigenvalues 
bundles in fig. ^. They can be numbered by q. We found {Kq, Xq) to 
be (2/73,3/2), (72/ v^, 73/2), ((^5 - l)/73, ^15 - 675/2), (2/3,1/2) 
for q = 3,4,5,6, respectively. Positions of the cross-points do not change 
with the lattice size. This structure of the cross-points will be discussed in 
more detail elsewhere. Here we only want to mention some features of this 
structure. All the points lie in 'antiferromagnetic' phase. It is seen in fig.^ 
that for g — i> oo the points approach the border to the 'ferromagnetic' phase 
Kq — i> Kq and that the corresponding eigenvalues Xq decrease. Moreover 
the probability of encountering an elementary loop of length q on random 
triangulation falls off exponentially @] which means that there are very few 
loops with high q. Taking all these facts into account one may expect that the 
presence of discrete line of cross-points which approaches Kq shall dominate 
the shape of the spectra for the lowest eigenvalues at Kq, leading in particular 
to the appearance of discrete peaks also in the limit N ^ oo. 

On the other hand, one also expects that the discreteness of the spectrum 
disappears when K goes deeper into the ferromagnetic phase. The crossover 
between the two regimes depends on some combination of AK = Kq — K and 
the system size N. For AK larger than zero, indeed the amplitude of the 
fluctuations decreases with the size as can be seen in fig.^. The fluctuations 
disappear in the limit N —^ oo leaving out a smooth distribution. 

The spectra Pj{x) of other eigenvalues Xj, j > 0, also exhibit the same 
fluctuating pattern at Kq as for j = 0, however the amplitude of the fluctua- 
tions decreases faster with AK and A^. As one can, for example, see in fig.0.c 
the oscillations are absent in the spectrum of the third lowest eigenvalue P2{x) 
already for K = 0.4727. 

For K far from 1/73 density distributions are smooth. In this case, the 

situation resembles the one known for instance from the considerations of 

chiral matrix models : separate densities Pj{x) are described by bell-shaped 

functions [^. They sum up p(x) = ^,Pj(a;) to a function with oscillations 

27| like for example in fig.|^.d. It would be interesting to find a random 



matrix model which reproduces ^j{x) and p(a;) analytically. 

Another interesting feature of the histograms is the presence of singular 
peaks for which the number of entries grow with the size of the lattice. The 
peaks lie outside the range displayed in the figures. We show in the appendix 
that such peaks are also present in the spectrum on the regular lattice. In this 
case we calculated analytically that the height of the peaks is logarithmically 
divergent in the lattice size. Thus this seems to be a generic situation. 
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Figure 6: The probability distribution po{x) for (a) Kq = l/VS , N = 64, 
(b) Ko = 1/VS , N = 128. (c) K = l/V^e-^-^ = 0.4727, A^ = 64 (d) 
K = l/y3e~°-2 = 0.4727, N = 128. The amphtude of the oscillations in the 
bottom row decreases with N. The histograms (a) and (b) contain 2 x 10^ 
counts each, (c) 6 x 10^, (d) 7 x 10^. The bin size is 5 x 10"^. 



Let us come back to the universal properties of the system of fermions 
interacting with gravity. We will study now in more detail behavior of the 
spectrum at the critical point of the Ising model. More precisely we shall 
be interested in the mass gap of the of the spectrum at K close to Kcr- We 
define the gap as the position of the center of mass for the distribution of 
the lowest eigenvalue po{x). We denote it by M = f xpo{x). We want to 
determine the dependence M = M{K) on the hopping parameter K for the 
given system size N. We do this numerically using the Lanczos algorithm^. 

^The Lanczos algorithm |Q is an iterative procedure to calculate eigenvalues. It is 
frequently used to approximately determine the lowest part of the eigenvalue spectra of 
large matrices, for which exact standard diagonalization algorithms would require a too 
long time. In a single iteration step the Lanczos algorithm finds one approximated eigen- 
value and improves quality of the previously calculated ones. As a rule it first produces 
the smallest and the largest eigenvalues and then successively fills up the remaining part 
of the spectrum, The accuracy increases with the number of iterations. We checked in 
our case using matrices of sizes up to 128, that when we keep the number n of iteration 
proportional to the size of the matrix n — cN with, c = 0.25, the distribution of the lowest 
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Figure 7: The figures show probabihty distributions Pj{x) of the lowest 
eigenvalues for A^ = 64 and (a) K = 1/V3 = 0.5774, (b) K = l/VSe~°-^ = 
0.5224, (c) K = l/y3e-°-2 = 0.4727, (d) K = l/v^e^°-3 = 0.4277. More 
precisely, each figure (a-c) contain distributions for j = 0,1,2, while the 
figure (d) for j = 0, 1, . . . , 9. Additionally, in the figure (d) the p.d.f. p{x) is 
shown. Each histogram Pj(x) contains 6 x 10^ counts and has the bin size is 
5 X 10~^, while p{x) contains 10^ counts, and has the bin size 1.25 x 10^^. 



For given N the function M{K) has a minimum (see fig.§). The value 
of the minimum M^, plays the role of a mass gap, while its position K^, of a 
pseudo-critical hopping parameter. We determined K^, and M^, for different 
system sizes. The results are collected in the table ||. We fitted the data 
points to the following finite size scaling formulas 



M. 



N-^H 



I 



t 

N 



K. 



Koo- 



N^ 



(29) 



The exponent d^ is the fractal dimension of the surface given. A typical 
linear extent of L scales as L = N^^'^" . When the physical mass is equal 
zero, L sets the scale for the correlation length. Its inverse gives the minimal 
eigenvalue of the spectrum M^,. For smaller systems one expects correc- 
tions to scaling. We took it into account by introducing a phenomenological 

eigenvalue agrees with the one obtained by an exact algorithm. 
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Figure 8: The position of the center of mass is shown for the distribution 
of the smallest eigenvalue of the operator T> for A^ = 96 as a function of the 
hopping parameter K. 



N 


K, 


M, 


32 

48 

64 

96 

128 

192 

256 

384 

512 

768 

1024 


0.352(3) 
0.360(2) 
0.364(2) 
0.368(2) 
0.370(1) 
0.372(2) 
0.372(2) 
0.374(2) 
0.374(2) 
0.375(3) 
0.375(1) 


0.1395(5) 
0.1162(4) 
0.1016(2) 
0.0858(2) 
0.0766(2) 
0.0656(3) 
0.0589(3) 
0.0509(3) 
0.0459(4) 
0.0397(2) 
0.0359(3) 



Table 1: Positions and values of the minima of the function M{K) repre- 
senting the center of mass of the distribution of the smallest eigenvalue of 
the operator V for different system sizes A^. 
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Figure 9: The mass gap M^ for different system sizes A^, and the curve 
representing the best fit to the formula (pOf) : l/dff = 0.348(4), b 
and t = 5.7(5). 



0.40(1) 



correction t/N to the formula (^). This correction significantly improves 
quality of the fit for the studied range of A^. The best fit to the formula is 
l/clff = 0.348(4), b = 0.40(1) and t = 5.7(5). The corresponding curve is 
plotted in fig.^. The curve fits indeed very well to all the data points. The 
error bars of the best fit parameters were estimated by jack-knife. 

We compared the goodness of the best fits to the formula (p9D and anal- 
ogous formulas in which the correction t/N was substituted by t/N^^"^ and 
t/N'^/"^. We obtained x^/d-o.f. = 0.52 for t/N while in the other two cases 
1.66 and 1.83, respectively. Thus, among those three correction types, the 
one t/N is best in this range. We have also checked that the fitted value 
l/du = 0.348(4) is stable against the successive removal of the data points 
of the smallest volumes. 

There are different theoretical predictions for the value of the Hausdorff 
dimension dn = 4: |23, 24 and (i// = 3 [25|, the latter of which was obtained 



for a test fermion in the gravitational background coupled to matter field 
with the central charge c = 1/2 ||25|. The Hausdorff dimension measured in 
our MC simulations dn = 2.87(3) favors dn = 3. One should, however, be 
aware that in measurements of the Hausdorff dimension there are large finite 
size effects, as can be seen from the considerations of pure gravity p3[ . 

The best fit for the second formula in (|29|) is given by Koo = 0.3756(16), 
K = 1.03(30) and a = 0.9(5) (see fig.pil|). The limiting value K^o is in agree- 
ment with the theoretically calculated critical value Kcr (0)- The scaling 
exponent k is almost equal 1 which would suggest a kinematic scaling saying 
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Figure 10: The pseudo-critical parameter K^ for different system sizes A^, 
and the curve representing the best fit to the scahng formula (p9D : K^o = 
0.3756(16), K = 1.03(30) a = -0.9(5), plotted as a function of 1/A^. 



that the average distance between eigenvalues decreases as 1/A^. 

Conclusions and outlook 

We have investigated the properties of the Dirac- Wilson operator on a ran- 
dom triangulation. In particular we have extracted from the spectrum of 
the operator the values of the fractal dimension dn and the critical value 
of the hopping parameter Kcr- At this value oi K = K^r the fermions be- 
come massless and one can take a continuum limit corresponding to massless 
fermions interacting with 2c? gravity. The value of K^r has been calculated 
analytically, however, du has not been unambiguously determined theoreti- 
cally. In the neighborhood of Kcr-, which lies deep in the ferromagnatic phase 
we found spectral distributions which are typical for random systems. 

Apart from Kcr there is another interesting value of K, namely Kq = 
l/Vo which lies on the border between the ferromagnetic and antiferromag- 
netic phases of the corresponding Ising model. We observed that the distri- 
butions of the lowest eigenvalues, Pj{x), becomes discrete when K goes to 
Kq from below. This is an unexpected phenomenon for a random system. 

For some values of the hopping parameter K in the 'antiferromagnetic' 
phase, eigenvalues of the Dirac- Wilson operator decouple from the random 
structure of the matrix and depend only on local geometrical properties of 
the triangulation. This, as we discussed, leads to the appearance of discrete 
spectra at Kq. 
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There are many natural extensions of the studies presented in this work. 
One should try to understand properties of the spectrum of the Dirac oper- 
ator from the point of view of the random matrix theory pB| , ^, ^ . This 
is a slightly different type of randomness than the one provided by the cou- 
pling to the vector gauge field which is usually discussed in the context of 
QCD. However exactly this type of randomness may be important in quan- 
tum gravity. 

Next, one can investigate in more detail the relation of the quenched ap- 
proximation to the full model. The quenched model describes a test particle 
in pure gravity. From this exercise one could perhaps draw a general lesson 
about the effect of quenching on the spectrum of the Dirac operator. This 
can be important because this type of approximation is frequently used in 
many physical contexts, for example, in QCD. However, one usually is not 
able to quantify the effects of quenching. 

Furthermore, one can study effects of changing topology by consider- 
ing non-spherical 2d-triangulations. As mentioned, this requires a careful 
treatment of various spin structures which may be admitted by a manifold. 
Contrary to higher dimensional case, where the existence of spin structure 



is related to the second Stiefel-Withney class |2^, here the question of the 
existence reduces to the orientability of the manifold. Also the classification 
of spin structures is relatively simple in the 2d case. The spin structures 
can be classified by a set of signs defined on all classes of non-contractable 
loops. The signs tell us whether boundary conditions for a parallel transport 
of a spinor around those loops are periodic (+1) or anti-periodic (—1). For 
a manifold with genus h, there are 2h different classes of non-contractable 
loops and hence there are 2"^^ different spin structures. 

Finally one should try to find a lattice implementation of the Dirac oper- 
ator for higher dimensional compact simplicial manifolds. Many parts of the 
construction can be directly generalized from the 2d case; actually almost 
all, except the link sign degrees of freedom, Sij (|^), which as it turns out are 
not sufficient in general case for the connections Uij to fulfill the consistency 
condition for all plaquettes (^). 
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Figure 11: Regular triangulation of the plane and its dual lattice. Fermions 
live on the vertices of the dual (honey-comb) lattice. The elementary cell 
contains two distinct dual node positions denoted by A and B. 



Appendix 

For a comparison in the appendix we calculate spectrum of the Dirac- Wilson 
operator on the regular planar triangulation built of equilateral triangles 
with fermion field located in the centers of triangles. If one connects the 
centers by links, they form a dual lattice; in this case it is a honey-comb 
lattice (see fig.^). It is convenient to divide the vertices of this lattice 
into two classes A and B forming a check-board. The fundamental cell 
on the triangulation contains one site of each. One reconstructs the entire 
triangulation translationally copying the fundamental cell using multiples of 
two the vectors di = Uq + ni, and d2 = Uq + n2 constructed from the link 
vectors no = (0, 1), rii = (^3/2, 1/2), na = (-^3/2, 1/2). 

Using translational symmetry of the lattice we can now rewrite the action 
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(0) in the following form 
K ^ - 

S = -—^^[(pi+dA'^ + nd-l)i^i,B + 'ipi,B{'i--nd--f)lpi+d,A\ 
i d=l 

K . 

- Y ^ [^i,A(l - no ■ -f)ipi,B + ^i,B(l + ^0 ■ 7)^*,^] 

i 

+ - ^ [i'i,Ai'i,A + i'i,Bi'i,B] , (30) 

i 

where the first index in ipi^A is a double index consisting of two integers 
{ii,i2), which give the position of the cell x = iidi + 12(^2, while the second 
label denotes the position A or B within the cell. In the component notation, 
the addition of di to i corresponds to [11^12) — ^ (^1 + 1,'^2), and of d2 to 
(^1,22) —^ {H,i2 + !)• In the expression ( ^0]) we have used a shorthand 
notation denoting the sum over di and (^2 by c/ = 1, 2. We can now partially 
diagonalize the problem using the Fourier transform of the index i = {ii, 12) 
to the momentum space p = (pi, ^2)- This leads us to a block-diagonal matrix 
consisting of four by four blocks. Each block D{p) corresponds to one Fourier 
mode ijjpD[p)iljp. The four by four matrix D{p) is indexed by the spinor index 
of ip and of the position label A or B. For each p, diagonalization of D{p) 
yields four eigenvalues 



_ 1 , ..v^ 

Ap — 

where 



Xp = -±K^^Jw±zy/A-{w^\ (31) 



w = cos(pi) + cos(p2) + cos(pi - P2) ■ (32) 

The distribution of eigenvalues (pTf) on a finite lattice L x L with periodic 



boundary condition in the (ii 2 directions is shown in fig.[T^. In this case the 
momenta admit the values pi^2 = Svrfci 2/i^, where ki^2 = 0, . . . , L — 1 and 
hence the operator has 4L^ eigenvalues. 

Similarly, one can find eigenvalues of the operator V 



Xp = ±^]J- + —iw + Q)±^i—iw -3) + 2^ + 9X^-4. (33) 
Using this formula we can calculate the spectral density 

p(x)=2im^$^5(a;-zA). (34) 
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Figure 12: Eigenvalues A of the Dirac- Wilson operator on a regular trian- 
gulation with L = 50 and K = 0.33. 



The spectrum terminates at a small positive value K^r (see fig.|TB|). It goes 
to zero only for K = K^ 



1 

3' 



In the large L-limit, the two peaks in fig. |13 
develop a logarithmic singularity. 

For the regular lattice the critical value of the hopping parameter is given 
by the standard equation K^r = 1/q', where q is the number of links emerging 
from the vertex. In this case g = 3. For random lattice this condition 
is dressed by lattice fluctuations. Although each vertex has coordination 
g = 3, the critical value of the hopping parameter is shifted from 1/3 to the 
value given by equation (p3[). For the regular lattice, the spectrum has an 
eigenvalue equal exactly zero for the critical value of the hopping parameter. 
This is not the case for random lattice, where the smallest eigenvalue has a 
distribution whose center of mass approaches zero only for large A^. On the 
regular lattice, the lowest part of the spectrum does not move when N goes to 
infinity, but becomes denser. The average distance between the eigenvalues 
scales like N~^^'^" , with the canonical dimension cij^ = 2 while on the random 
lattice, the position of the lowest eigenvalue moves towards zero with N~^''^" 
with a dressed exponent du = 2.87(3) resulting from the fractal structure of 
fluctuating geometry. 
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Figure 13: Histograms of eigenvalues of the operator T> for L = 3000 and 
K = 0.3. At the critical value K = K^r the spectrum goes continuously to 
zero, while for K ^ Kcr, like for example for K = 0.3 presented in the figure, 
its low- A part is cut off at some positive Xmin- 
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